!> 测试版本

module exp

    use gauss, only: swap

contains

    !> code:blocks群友提供。
    pure subroutine solve(a, b)
        ! n*n * n*m = n*n
        real(kind=8), intent(in) :: a(:, :)
        real(kind=8), intent(inout) :: b(:, :)

        real(kind=8) :: C(size(A, 1), size(A, 2) + size(b, 2))
        integer      :: n, m
        real(kind=8) :: temp, temp2
        integer      :: i, j, k, major

        n = size(a, 1)
        m = size(b, 2)
        C(:, :n  ) = a
        C(:, n+1:) = b

        ! column major method
        do i = 1, n - 1
            ! 寻找第i列主元 find major position
            major = maxloc(abs(c(i:n, i)), dim=1) + i - 1
            !数据交换,cahce miss
            call swap(c(i, i:), c(major, i:))
            if (abs(c(i, i)) < epsilon(temp)) return
            temp = 1.d0/c(i, i)
            !计算
            do j = i + 1, n+m
                c(i, j)       = c(i, j)*temp!归一化
                c(i + 1:n, j) = c(i + 1:n, j) - c(i + 1:n, i)*c(i, j)
            end do
            !第i列
            c(i, i)       = 1.d0
            c(i + 1:n, i) = 0.d0
        end do

        temp    = 1.d0/c(n, n)
        c(n, n) = 1.d0
        ! 回代
        do j = 1, m
            !归一化
            b(n, j) = b(n, j)*temp
            do k = 1, n - 1
                temp2         = b(n - k + 1, j)
                b(1:n - k, j) = b(1:n - k  , j) - temp2*c(1:n - k, n - k + 1)
            end do
        end do

    end subroutine solve

    !> code:blocks群友提供。
    pure function solve1(a, b) result(x)
        ! n*n * n*m = n*n
        real(kind=8), intent(in) :: a(:, :)
        real(kind=8), intent(in) :: b(:, :)
        real(kind=8) :: x(size(b, 1), size(b, 2))

        real(kind=8) :: C(size(A, 1), size(A, 2) + size(b, 2))
        integer      :: n, m
        real(kind=8) :: temp, temp2
        integer      :: i, j, k, major

        n = size(a, 1)
        m = size(b, 2)
        C(:, :n  ) = a
        C(:, n+1:) = b

        ! column major method
        do i = 1, n - 1
            ! 寻找第i列主元 find major position
            major = maxloc(abs(c(i:n, i)), dim=1) + i - 1
            !数据交换,cahce miss
            call swap(c(i, i:), c(major, i:))
            if (abs(c(i, i)) < epsilon(temp)) return
            temp = 1.d0/c(i, i)
            !计算
            do j = i + 1, n+m
                c(i, j)       = c(i, j)*temp!归一化
                c(i + 1:n, j) = c(i + 1:n, j) - c(i + 1:n, i)*c(i, j)
            end do
            !第i列
            c(i, i)       = 1.d0
            c(i + 1:n, i) = 0.d0
        end do

        temp    = 1.d0/c(n, n)
        c(n, n) = 1.d0
        ! 回代
        do j = n+1, n+m
            !归一化
            c(n, j) = c(n, j)*temp
            do k = 1, n - 1
                temp2         = c(n - k + 1, j)
                c(1:n - k, j) = c(1:n - k  , j) - temp2*c(1:n - k, n - k + 1)
            end do
        end do

        x = c(:, n+1:)

    end function solve1

end module exp